home *** CD-ROM | disk | FTP | other *** search
/ FM Towns: Free Software Collection 7 / FM Towns Free Software Collection 7.iso / t_os / pigcc / pigcc.doc < prev    next >
Text File  |  1993-11-30  |  13KB  |  248 lines
  1. ＧＣＣ＆ＧＡＳによる円周率の計算　　　　　　　　　　　Ｂｙ　Ｊｏｕｊｉ
  2.  
  3. ★まえがき
  4. 　私が、FM ＴＯＷＮＳ ＵＸ２０を買ったのは１９９２年２月初めでした。そして２
  5. 月の末に、迷ったあげくＧＮＵ ＣＤ ｒｅｌ．２を注文しました。ところが、品不足
  6. という事でなかなか届かず、半分諦めかけていたのですが、５月の末になってようや
  7. く、待望のＧＮＵがやってきました。早速、ＧＣＣでお絵描きプログラムを作ったり、
  8. あっちこっちのサンプルプログラムを再コンパイルしたりしてみました。やっぱり、
  9. 自分で作ったプログラムがＴＯＷＮＳで走っているのを見るのは最高の気分ですね。
  10. 　そのうちに、Ｏｈ！ＦＭ１９９０年７月号の若松登志樹氏の記事にも出ていた円周
  11. 率の計算をしてみようと思い立ちました。まずは、フリコレ３に入っていた片山慎太
  12. 郎氏作のＢＡＳＩＣプログラムを見て計算方法を勉強し、ついでに計算をしてみたと
  13. ころ、１００桁を２８秒で計算できました。そのプログラムを多少自分なりに改良し
  14. てみたら１００桁で２０秒になりました。さらにＣに移植してＧＣＣでコンパイルし
  15. たら、１００桁を０．５秒で計算できました。さすがＣは速いなと感心し、その時か
  16. ら、どの位まで速くできるかという挑戦が始まりました。
  17.  
  18.  
  19. ★円周率計算プログラム
  20. 　最初はＣで作っていたのですが、Ｃでは６４ビットの整数が無いため、３２ビット
  21. ×３２ビットおよび６４ビット÷３２ビットの演算が簡単ではありません。やはり、
  22. このような演算ルーチンはアセンブラで作らなければなりません。３８６ＣＰＵはこ
  23. のようなビット数の演算は１命令で実行できるので、アセンブラで作ればかなりの高
  24. 速化が期待できます。ＧＮＵのアセンブラはＧＡＳですが、ＧＡＳにはマクロ定義機
  25. 能などは無いようです。そこで、ＧＡＳのソースプログラムをプリプロセッサＣＰＰ
  26. に通して、その後アセンブルするようにしました。ＣＰＰの置換機能を利用すればＧ
  27. ＡＳもなかなか使えるものです。具体的には、 make.bat により、
  28.  
  29. 　  make  pigcc.s  2sen  exec
  30.  
  31. のようにコマンドラインから入力すると、ＣＰＰにより pigcc.s から pi2sen.s を
  32. 作成し、その pi2sen.s に対してアセンブル、リンク、実行（．ＥＸＰ）ファイルの
  33. 作成を行い、続いて実行します。最後の exec パラメータが無ければ実行はしません。
  34. なお、このバッチの実行には環境変数ＧＮＵが必要です。 pigcc.s は、もともとは
  35. Ｃプログラムをコンパイルしてアセンブラソースを出力させたものに手を加えていっ
  36. たものですが、今となってはほとんど原型をとどめていません。ｍａｉｎ以外は完全
  37. にオリジナルです。ライブラリ関数はｐｒｉｎｔｆとｃｌｏｃｋを使用しています。
  38. 　pigcc.s は２０００桁計算用になっていますが、ＮＤＷの定義値を１行目のコメン
  39. トに従って変更することにより計算桁数を変更できます。コメント文では計算誤差を
  40. 見込んで１０桁余分に設定しています。また、定義値ＦＤＳは、出力制御用のスイッ
  41. チであり、０にすると計算した全桁を出力しますが、１にすると最初の７２桁だけ出
  42. 力します。これは、計算ルーチンを高速化している過程で使用していたもので、計算
  43. 時間だけ知りたい場合に有効です。何しろ１万桁も計算すると画面出力の時間だけで
  44. もかなりのものですから。
  45. 　プログラムの実行は、Ｔ－ＯＳ Ｖ２．１Ｌ１０，２０のコンソール環境で次のよ
  46. うに行います。
  47.  
  48. 　　run386  pi2sen.exp
  49.  
  50. この場合、計算結果を画面に出力します。計算結果は改行も空白も挿入せず、べたに
  51. 出力されます。そして、最後に計算時間を出力します。start は計測開始時間、c_en
  52. d は２進数での計算終了時間、d_end は計算結果の出力終了時間です。時間の単位は
  53. １／１００秒です。
  54. 　出力をファイルにリダイレクトする場合は、後で述べる ketaori.exe を使って、
  55.  
  56. 　　run386  pi2sen.exp | ketaori >pi2sen.dat
  57.  
  58. というように、桁折りしながら出力した方が良いと思います。
  59. 　さて、このプログラムのアルゴリズムですが、マチンの公式と呼ばれている次式
  60. （１）を使っています。
  61.  
  62. 　　π　　　　　　　　 １　　　　　　　 １
  63. 　　―　＝　４・arctan ―　－　arctan ―――　　　　　　　‥‥‥（１）
  64. 　　４　　　　　　　　 ５　　　　　　 ２３９
  65.  
  66. アークタンジェントarctan（ｘ）は、グレゴリー級数と呼ばれている次式（２）の様
  67. なべき級数に展開できます。
  68.  
  69. 　　　　　　　　　　　　　　１　　　　　　１
  70. 　　arctan（ｘ）＝　ｘ　－　―ｘ＾３　＋　―ｘ＾５　・・・
  71. 　　　　　　　　　　　　　　３　　　　　　５　　　　　　　‥‥‥（２）
  72.  
  73. 式（１）の arctan を式（２）のべき級数に置き換えることによって、πが計算でき
  74. ます。Ｏｈ！ＦＭの若松氏の記事ではシュテルマーの公式を使用していましたが、私
  75. が計算してみた結果ではマチンの公式の方が少し速かったようです。その差は僅かで
  76. したが。
  77.  
  78.  
  79. ★計算速度について
  80. 　計算時間は、目標としていた若松氏のものより高速になり、一応、目標達成といっ
  81. たところです。本プログラムと若松氏のプログラムの計算時間の比較を【表１】に示
  82. します。若松氏のプログラムはフリコレ５に入っていたものです。出力は全てＲＡＭ
  83. ディスク上のファイルにリダイレクトしています。
  84.  
  85. 【表１】　　円周率の計算時間（FM ＴＯＷＮＳ）
  86. +----------+----------------+----------------+----------------+
  87. |          | 本プログラムを | 若松氏のものを | 若松氏のものを |
  88. |          | ＵＸ２０で実行 | ＵＸ２０で実行 | モデル２で実行 |
  89. +----------+----------------+----------------+----------------+
  90. |   １万桁 |         ５２秒 |     １分２１秒 |     １分１５秒 |
  91. |          |       (48 + 4) |      (53 + 28) |                |
  92. +----------+----------------+----------------+----------------+
  93. | １０万桁 |   １時間２５分 |   ２時間０２分 |   １時間５２分 |
  94. |          |   (4776 + 325) |  (5310 + 2002) |                |
  95. +----------+----------------+----------------+----------------+
  96.  
  97. 　【表１】の「モデル２で実行」というのは、若松氏の結果をそのまま載せたもので
  98. す。これを見ると、このような数値計算はやはり３８６ＤＸが速いのですね。１６Ｍ
  99. Ｈｚノーウェイトにもかかわらず、３８６ＳＸでは１０％弱遅くなっています。
  100. 　【表１】の下段のかっこ内の数値は、最初のものが２進数での計算時間、次のもの
  101. が２進→１０進変換および出力時間（単位は秒）です。本プログラムでは、２進→１
  102. ０進変換および出力時間がかなり短縮されているのが分かると思います。若松氏のプ
  103. ログラムでは円周率の２進計算結果に対して１０倍しては整数部分を出力するという
  104. ループになっています。本プログラムでは、２進計算結果に対して１０００００００
  105. ００倍（３２ビットで表せる最大の１０のべき乗数＝１０＾９）しては整数部分を２
  106. 進→１０進変換して出力するというループになっています。これによって、円周率の
  107. 多桁数値に対するかけ算の回数が１／９になります。
  108. 　また、本プログラムでは、２進→１０進変換して出力するという部分は、実際には
  109. Ｃライブラリのｐｒｉｎｔｆ関数の％ｄ変換を行っているにすぎません。この部分も
  110. 自前のルーチンを用意したらさらに高速になるのではないかと思われる方もいるでし
  111. ょう。私もそう思い自前の変換ルーチンをつくり、ＤＯＳのファンクションコール９
  112. 番によって９文字ずつ出力するようにしてみました。こうすると、画面出力する場合
  113. はｐｒｉｔｆよりも高速になりますが、出力をＲＡＭディスク上のファイルにリダイ
  114. レクトした場合はなんとｐｒｉｎｔｆの方が高速であることが分かりました。
  115. 　この理由ですが、ｐｒｉｎｔｆ関数は多分自前のバッファを持っていて、リダイレ
  116. クトした場合それが有効に機能するのではないかと思います。ＤＯＳのシステムコー
  117. ルを呼んだときに生じるｒｕｎ３８６．ｅｘｅのプロテクトモード←→リアルモード
  118. 自動変更機能は、便利ではありますがかなりのオーバーヘッドを生ずると考えられま
  119. すから、自前のバッファを持ちシステムコールの回数を減らせばかなりの高速化が期
  120. 待できます。ｐｒｉｎｔｆがそれをやってくれているのならば、簡単だし高速だしの
  121. 一石二鳥であり、わざわざ自前の出力ルーチンをつくる必要は無いわけです。
  122.  
  123.  
  124. ★ＧＡＳの表記法について
  125. 　ここで、簡単にＧＡＳの表記法について説明します。ＧＡＳの表記は、ＭＡＳＭな
  126. どとはかなり異なっている点がありますから注意して下さい。いちばん混乱し易いの
  127. は、オペランドの順番がソース、デスティネーションの順であり、ＭＡＳＭなどとは
  128. 逆になっている点です。以下に、相違点をまとめます。
  129. 　●ニーモニックの末尾に、オペランドのデータサイズを示す文字ｂ／ｗ／ｌを
  130. 　　付ける。ｂ／ｗ／ｌは、バイト／ワード／ダブルワードを表わす。
  131. 　●オペランドの順番は、ソース、デスティネーションの順である。
  132. 　●レジスタの名前には、先頭に％を付ける。
  133. 　●即値オペランドには、先頭に＄を付ける。
  134. 例えば、
  135.  
  136. 　　movl %eax,%ebx
  137.  
  138. は、３２ビットのＥＡＸレジスタからＥＢＸレジスタへの転送を表わします。
  139.  
  140.  
  141. ★３８６プログラムの高速化について
  142. 　では、pigcc.s の中で高速化に留意した点を述べてみましょう。pigcc.s の中で、
  143. 計算時間に最も影響を与えるのが、サブルーチン＿ｄｉｖ＿ａｄｄと＿ｄｉｖ＿ｓｕ
  144. ｂです。この２つのサブルーチンはほとんど同じ構造になっていますから、＿ｄｉｖ
  145. ＿ａｄｄによって説明しましょう。＿ｄｉｖ＿ａｄｄの中でも、
  146.  
  147. 　　L40:
  148. 　　          ･
  149. 　　          ･
  150. 　　        loop L40
  151.  
  152. のループ内ルーチンが実行時間に最も影響するので、この中を重点的に高速化しまし
  153. た。７つの汎用レジスタは全て変数または定数に割り当てて除算２回と加算１回を連
  154. 続的に実行しています。工夫した点は、ｘｃｈｇ命令を使ってレジスタの内容を交換
  155. し、１ループ内で２回の除算を行うようにしたことです。
  156. 　また、このループ内の
  157.  
  158. 　　        adcl $0,-4*NDW-8(%esi)
  159.  
  160. にも注目してください。この行はなくても、同じ動作をするようにプログラムする事
  161. が出来ますが、この行を入れておくことにより条件分岐せずに処理が進む確率が圧倒
  162. 的に高くなります。３８６はプリフェッチ、プリデコードをパイプラインで処理して
  163. いますが、分岐命令があるとプリフェッチ、プリデコードの結果が無効になってしま
  164. い、処理時間が遅くなってしまいます。できるだけ条件分岐や、ｊｍｐ命令を使わな
  165. いようなプログラムにする事が大切です。
  166. 　そして、すぐ次の行の
  167.  
  168. 　　        jc L42
  169.  
  170. ですが、この飛ばし方と飛ばした先からの戻し方もちょっと変ですね。今説明したよ
  171. うに、すぐ上の行のおかげで、この行の条件分岐が成立する確率は非常に小さくなっ
  172. ています。条件分岐命令は、分岐しないときの方が速いですから、このように命令を
  173. 配置することによって、条件分岐することなく高速に処理を進めることが出来ます。
  174. また、条件分岐する確率は非常に小さい（約４０億分の１）のですから、分岐先のル
  175. ーチンはそれほど高速化に気を使わなくても良いわけです。
  176.  
  177.  
  178. ★高速化のまとめ
  179. 　私がこのプログラムを作成し高速化するに当たって気を付けた点を、以下にまとめ
  180. てみます。
  181. （１）データ処理単位を大きくする。８ビットより１６ビット、１６ビットより
  182. 　　３２ビットにする。
  183. （２）複数のループで処理しているものを、１つのループにまとめられないかを
  184. 　　考える。
  185. （３）繰り返し回数の多いループ（内側のループ）を集中的に高速化する。
  186. （４）データのロード，ストア，転送の回数を減らす。
  187. （５）データは可能なかぎりレジスタに保存する。レジスタに対する変数の割り
  188. 　　付けを最適にする。
  189. （６）レジスタに割り付けることが出来ない変数は、高速なアドレッシングモー
  190. 　　ドでアクセスできる格納場所を選ぶ。
  191. （７）分岐を少なくする。条件分岐を入れる場合は、頻度の少ない場合に分岐す
  192. 　　るようにする。
  193.  
  194.  
  195. ★その他のの注意点
  196. 　ＧＣＣでコンパイルしたアセンブラソースを見ていて気が付いたのですが、ＧＣＣ
  197. コンパイラは、レジスタｅａｘ，ｅｃｘ，ｅｄｘをテンポラリな変数として割り当て
  198. るようです。すなわち、これらのレジスタは関数やサブルーチンの内部でその値を保
  199. 存する必要がないわけです。それ以外のレジスタは値が保存される通常の変数として
  200. 割り当てられるようですから、関数やサブルーチン内で値を変更する場合は、その入
  201. 口で保存し、出る前に復帰しておかなければなりません。また、ｅｂｐはスタックフ
  202. レームを指すポインタとして使われていますから、あまりいじらない方が良いでしょ
  203. う。
  204. 　また、これは円周率の計算時間を測定しているときに気が付いたのですが、Ｔ－Ｏ
  205. ＳでＯＡＫを組み込んだときと、ＯＡＫを外したときでは、計算時間が変化するので
  206. す。最初は原因が分からず、計算時間がころころ変わるので不思議だったのですが、
  207. ＯＡＫを組み込んだり外したりしているせいだと分かりました。ＯＡＫを組み込むと、
  208. １％程計算時間が長くなります。多分、ＯＡＫの割り込み処理か何かが原因だろうと
  209. 思います。ということは、ＯＡＫを外すとＴＯＷＮＳが１％速くなるということでし
  210. ょうか。なお、表１の計算時間の測定はＯＡＫを外して行っています。
  211.  
  212.  
  213. ★桁折りプログラム
  214. 　リスト１のプログラムは計算した全桁数を改行も何もしないで、べたうちに出力し
  215. てしまいます。私は、エディタはＶｚを使っていますが、このプログラムで円周率１
  216. ０万桁を計算し、その出力ファイルをＶｚで確かめようとしました。そうしたら、ど
  217. うしたことでしょう、突然Ｖｚがハングアップしてしまいました。考えてみると、１
  218. 行が１００ＫＢｙｔｅにもなっているのです。ということは６４ＫＢｙｔｅの壁を越
  219. えており、ハングアップしても無理はないなと思いました。
  220. 　ということで、長大な行を任意の桁数で桁折りするプログラムも作りました。keta
  221. ori.c がそれです。何の変哲もないフィルタプログラムで、標準入力と標準出力を使
  222. います。パラメータで桁折りする桁数を指定します。パラメータがなければデフォル
  223. トの１００桁で桁折りします。私は、ＬＳＩ Ｃ－８６試食版でコンパイルしてＥＸ
  224. Ｅファイル ketaori.exe にしています。次のようにパイプで使用すると良いでしょ
  225. う。
  226.  
  227. 　　run386  pi10man.exp | ketaori >pi10man.dat
  228.  
  229.  
  230. ★あとがき
  231. 　私は、この円周率計算プログラムで１０万桁までは計算して確認していますが、１
  232. ００万桁の計算は実行していません。このプログラムでは計算時間は桁数の２乗に比
  233. 例しますから、１００万桁の計算は５．９日もかかります。計算途中経過をセーブで
  234. きるようにして、細切れに計算できるようにすれば良いと思いますが、なかなか面倒
  235. なものですから。
  236. 　円周率計算プログラムも作り終わった今頃になって、大野栄一氏の「パソコンで挑
  237. む円周率」と金田康正氏の「πのはなし」を読みました。やっぱり円周率って奥が深
  238. いですね。また機会があったら、ガウス・ルジャンドルの公式にも挑戦してみたいと
  239. 思います。
  240. 　最後になりましたが、この記事が、読者諸兄がプログラム開発するときに、少しで
  241. もお役に立てば幸いです。
  242.  
  243. ★参考文献
  244. 若松登志樹「ＹＥＳ，Ｉ ＨＡＶＥ Ａ ＮＵＭＢＥＲ」Ｏｈ！ＦＭ１９９０年７月号
  245. 安藤寿茂（ｅｌｆｉｎ）「ｇａｓマニュアル」ＧＮＵ for ＴＯＷＮＳ Rel.2
  246. 大野栄一「パソコンで挑む円周率」講談社，ブル－バックス
  247. 金田康正「πのはなし」東京図書
  248.